Adaptive seismic noise and interference attenuation method

ABSTRACT

A method relating to filtering coherent noise and interference from seismic data by constrained adaptive beamforming is described using a constraint design methodology which allows the imposition of an arbitrary predesigned quiescent response on the beamformer. The method also makes sure that the beamformer response in selected regions of the frequency-wavenumber space is entirely controlled by this quiescent response, hence ensuring signal preservation and robustness to perturbations. Built-in regularization brings an additional degree of robustness. Seismic signals with arbitrary spectral content in the frequency-wavenumber domain are preserved, while coherent noise and interference that is temporarily and spatially nonstationary is adaptively filtered. The approach is applicable to attenuation of all types of coherent noise in seismic data including swell-noise, bulge-wave noise, ground-roll, air wave, seismic vessel and rig interference, etc. It is applicable to both linear and a real arrays.

This invention relates to seismic data acquisition and to methods ofprocessing seismic data. It relates to a process for filtering coherentnoise and interference from seismic data by an adaptive beamformingmethod. In another aspect, it relates to adaptively filtering coherentnoise and interference from seismic data while preserving seismicsignals with arbitrary spectral content in the frequency-wavenumberdomain. In yet another aspect, it relates to adaptively filteringcoherent noise and interference that is temporally and spatiallynonstationary. In a further aspect, it relates to adaptively filteringcoherent noise and interference that has been recorded by a sensor arrayin the presence of perturbations.

BACKGROUND OF THE INVENTION

In seismic surveys, a seismic source induces seismic waves at or nearthe surface of the earth. Explosives, vibrating devices and airguns areexamples of seismic sources. These waves propagate through the earth andare reflected, refracted, and diffracted by formations within the earth,and can be detected by a plurality of sensors (typically geophones orhydrophones) at the earth's surface. Each such receiver monitors theseismic wavefield, which is then recorded. The data received by areceiver and then recorded are collectively called a trace. Thecollection of traces are stored for further processing to gaininformation about the earth's subsurface. Such information is commonlyinterpreted to detect the possible presence of hydrocarbons, or tomonitor changes in hydrocarbon bearing rocks.

Seismic data in general contains coherent noise signals, along withseismic reflection signals. These noise signals, hereafter referred toas the noise, interfere with the interpretation of the seismic signals,and degrade the quality of the subsurface images that can be obtained byfurther processing. It is therefore desirable to suppress the noise thatis present in the recorded data before processing it for imaging.

In land seismics, source-generated noise like ground-roll and air-wavesare the dominant noise types, and can lead to severe degradation in dataquality. In marine seismics, energy propagating as waves trapped in thewater column and near-surface layers is a significant source, as well asswell noise and bulge-wave noise which result from waves propagatingalong the streamers of receiver devices. Other sources of coherent noisein marine seismics include passing vessels, other vessels acquiringseismic data in the vicinity, or any drilling activity close to thesurvey area.

An important feature of the so-called coherent noise present in seismicdata is the distance over which the noise appears coherent. In manycircumstances, the noise is coherent over only a few meters. In othercases, although the noise is mostly coherent, there exists spatiallyimpulse noise. in such cases, filtering methods which have large spatialextent, like the known frequency-wavenumber filtering generateundesirable artifacts, which are mistakenly identified as seismic eventsafter further processing and imaging.

Another feature of the noise present in seismic data is that it is oftennon-stationary as a function of time; i.e. its characteristics change asa function of time.

During recent years there have been suggested a variety of methodsemploying the central concept of applying adaptive signal processingideas to the problem of suppressing coherent noise in seismic data.Booker and Ong, in: “Multiple-constraint adaptive filtering,”Geophysics, Vol.36, pp. 498-509, 1971, have derived an algorithm formulti-channel time-series data processing, which maintains specifiedinitial multiple filter constraints for known signal or noise sourceswhile simultaneously adapting the filter to minimize the effect of theunknown source field. The constraints are of the multiple look-directionconstraints type, where look-directions must be precisely specified.

In the International Patent Application WO97/25632, Ozbek has describeda class of adaptive signal processing techniques for attenuation ofdispersive, nonstationary and aliased coherent noise in seismic data, inthe presence of phase and amplitude perturbations. The methods developedcan be classified as multi-channel adaptive interference cancellers.Since a signal-free noise reference is not readily available in seismicdata acquisition, various preprocessing techniques are introduced togenerate the coherent noise reference channels. In the single-componentversion of the method, moveout (apparent velocity) and spatio-temporalcoherence are used as the criteria for differentiating between thesignal and the noise. In the multi-component version, polarization isused as an additional attribute for differentiation. Once single ormultiple noise reference channels are established, coherent noise in theprimary channel is canceled using data-adaptive least-squaresmulti-channel filter banks.

U.S. Pat. No. 4,556,962 attempts to attenuate the ground roll from asurface seismic source by placing a sensor close to the source to detectthe interfering noise. The interfering noise is scaled, delayed andsummed with signals from a more distant geophone array and thencross-correlated with the original vibrational source. This patent alsosuggests that an adaptive filter be used so as to modify the delayedsignal to correspond more closely to that detected by the more distantgeophone array. However, ground roll is in general of an inhomogeneousnature; due to dispersion and scattering from near surface anomalies theground roll measured at one point increasingly deviates in characterfrom that measured at another with increasing distance. Hence, theground roll measured close to the source may be substantially differentfrom that received by the geophone array, and the adaptive filter maynot be able to deal with this. It is also difficult to measure seismicsignals (ground roll) close to the source. Often the nearest offset is100 meters. For close measurements, more robust sensors may be neededand detector ‘character’ matching should be an important preliminarystep.

In U.S. Pat. No. 4,890,264 a method is proposed for suppressingnon-uniformly distributed noise generated by surface wave propagation,wind, and machinery. A number of horizontally sensitive geophones aredistributed amongst the conventional vertically oriented geophones. Theoutputs of the surface wave detectors are utilized in conjunction withan adaptive filter to cancel the effects of the surface waveinterference. This method for the suppression of ground roll isinherently a multi-component method, and cannot be used in conjunctionwith single component acquisition. In addition, it neglects the factthat some seismic body wave energy also is detected by horizontallysensitive geophones, and this may cause signal cancellation.

In UK Patent Application GB-A-2273358 linearly constrained adaptivebeamforming and adaptive interference canceling beamforming is used forground roll suppression. In linearly constrained adaptive beamforming,signals measured by an array of geophones are filtered and summed so asto preserve signals incident from a preferred direction whilesuppressing interferences incident from other directions. In applyingadaptive interference canceling, the moveout differential between theseismic reflections and the ground roll is used to form primary andreference channels. The filtering is performed using a continuouslyadaptive method such as the LMS (least-mean-square) algorithm. Thesuggested application is in seismic while drilling, where the horizontaloffset range is very small, so that the seismic reflections have analmost vertical angle of incidence and there is effectively a lot ofdata available from each source-receiver position since the roller conedrill bit used as the seismic source moves very slowly. The statisticsof the noise then change very slowly, allowing stochastic gradient typeof algorithms like the LMS to converge. However, in surface seismicexperiments the ground roll present is often non-stationary andinhomogeneous. Therefore, stochastic gradient type of algorithms such asLMS may be too slow in converging.

U.S. Pat. No. 5,237,538 proposes a method for removing coherent noisefrom seismic data. In this method, first the moveout characteristics ofthe noise are identified. Next, a space-time gate containing the noisedefined and extracted, and the moveout removed to flatten the noisetrain. Amplitudes and time variations are then removed from the gate.The coherent noise is estimated using conventional stacking. Asingle-channel Wiener filter is used to match the noise estimate to thenoise in the data trace containing signal-plus-noise. Having subtractedthe filtered noise estimate, inverse amplitude scalars are applied toundo the effect of amplitude equalization. The signal is then moveoutrestored into the original seismic record. This particular method forremoving coherent noise from seismic data is an application of thewell-known technique called Postbeamformer Interference Cancelling. Ithas some particular shortcomings for application for ground rollattenuation. First, the signal always leaks into the ground rollestimate, especially for shorter arrays. There is always a component ofthe signal present at the reference channel which is co-located in timewith the signal in the primary channel. On the other hand, when thearrays are allowed to be longer, the dispersion present in the groundroll make it difficult to achieve effective beamsteering.

In marine seismic surveys, an acoustic source generates waves whichtravel through the water and into the earth. These are then reflected orrefracted by the sub-surface geological formations, travel back throughthe water and are recorded by long hydrophone arrays which are towednear the surface of the water behind a seismic vessel. The hydrophonesare mounted in streamer cables, or streamers. There are usually 1-12streamers towed which are each several kilometers long. The streamersare made up of sections which may typically be 100-200 meters long; eachsection consists of hydrophones inside an outer skin which may be filledwith oil, foam, or a more solid substance. Stress-wires and spacers formthe internal skeleton of the streamer.

While the streamers are being towed behind the vessel, self-noise isgenerated due to a variety of sources. The lurching of the vessel,especially in rough seas, causes vibrations in the stress-wires whichinteract with the connectors and the oil-filled skin, generating bulgewaves (or breathing waves) which propagate down the streamers. Thepressure variations are detected by the hydrophones, adding andcorrupting the detected seismic signals. As the streamer moves throughthe water, boundary layer turbulence causes pressure fluctuations at theouter skin wall, which are again coupled to the hydrophones.

Bulge waves may also be caused by eddy shedding under elliptical watermotion about the streamer caused by wave action. A method of applyingadaptive signal processing to the attenuation of bulge waves isdescribed U.S. Pat. No. 4,821,241. There it is proposed to co-locatestress sensors with the hydrophones in the streamer. The stress sensorsare responsive to mechanical stresses applied to the cable, but aresubstantially unresponsive to acoustic waves propagating in fluid media.The signal outputs from the stress sensors are combined with the signaloutputs from the corresponding co-located hydrophones to cancel spurioussignals due to bulge waves.

Another method of applying adaptive signal processing to the attenuationof bulge waves was described is described U.S. Pat. No. 5,251,183. Inthis patent it is proposed to use an accelerometer secured between thelead-in section of the streamer and the hydrophone. Intra-shot andinter-shot accelerometer and hydrophone signals are recorded. The methodutilizes inter-shot and intra-shot adaptive processing loops. Theinter-shot adaptive processing loop derives inter-shot complex weightsfrom inter-shot accelerometer signals and inter-shot hydrophone signals.The intra-shot adaptive processing loop models bulge wave noise in theintra-shot hydrophone signals by combining the inter-shot complexweights with intra-shot accelerometer signals. Bulge wave noiseattenuation is achieved by subtracting the intra-shot bulge wave noisemodel from the intra-shot seismic detector signals.

SUMMARY OF THE INVENTION

In accordance with the present invention, there is provided a method forfiltering noise from discrete noisy seismic signals, comprising thesteps of receiving signals using a plurality of receivers; determiningpropagation characteristics of the signals with respect to receiverlocations; and filtering received signals using an at least partiallyadaptive filter such that signals having propagation characteristicsother than the determined propagation characteristics are attenuated.The filtering step comprising the step of defining at least twoindependent sets of conditions (constraints) with a first set defining adesired (quiescent) response and a second set defining the propagationcharacteristics of signals to be preserved and the step of adaptingfilter coefficients of the filter subject to the independent sets ofconditions (constraints) so as to minimize (optimize) the filter outputfor signals with propagation characteristics other than the determinedpropagation characteristics.

For the application of the invention, it is advantageous to define forthe optimization process of the filter weights or coefficients asignal-dependent part (correlation matrix) and a signal-independentpart. The signal-independent part usually comprises the constraints andis there often referred to as constraint matrix. Using this concept of aconstraint matrix, an important aspect of the invention can be describedas having within the constraint matrix a subspace which is defined bythe desired quiescent response and one subspace which defines theregions of the protected signal. By making these two subspacesorthogonal, filter weights can be found which simultaneously impose bothrestrictions upon the filter response. As the constraint matrixeffectively reduces the degrees of freedom of the filter available forthe adaptation process, this aspect of the invention can be described assplitting the total number of degrees of freedom into a first part whichis available for the adaptation process and a second part which is usedto define the constraints. The degrees of freedom assigned to theconstraints are split among those constraints which defines the desiredresponse and a second set defining the temporal and/or spatial spectralcontent or the propagation characteristics of the signals to bepreserved.

It is an advantage of the method to be not confined to narrow-bandsignals, but also applicable to wide-band seismic signals resulting in afilter that changes its response with the frequency of the inputsignal(s).

It is an important aspect of the invention that having derived a methodof separating the desired quiescent response and the constraintsrelating to the region into orthogonal subspaces, any known method canbe used for the adaptation process. Such known adaptive methods areknown and described in the literature, e.g. LMS, RLS, LSL, FTF, etc.

According to a preferred embodiment of the invention, a filter bankcomprising temporally and spatially local filters is used as theadaptive filter.

A filter bank can be defined as comprising M local multichannel adaptivefilters with K channels, each of a length L. For most applications, thenumber L of coefficients is equal to or larger than three. The number ofchannels K and of individual filters M are preferably two or more.

The use of a filter bank for noise attenuation of seismic signals hasbeen described in the International Patent Application WO97/25632.However, the present invention does not require defining a referencechannel in order to calculate the adapted filter bank coefficients. Nonoise estimate enters the adaptation process. Therefore, the presentmethod can be applied to noise contaminated seismic signals, where thereis no independent measurement or estimate of the noise available.

According to one aspect of the invention, the coefficients of the filterare constrained such that its response corresponds to that of abeamformer with a specified look-direction.

According to another aspect of the invention the constraints are setsuch that the filter preserves signals from a range of look directionsor of defined regions of the frequency-wavenumber domain. The region canbe pre-selected depending on the nature, more specifically on theapparent velocity of the seismic signals. Certain limits of thevelocity, such as 1500 m/s, define regions in the frequency-wavenumberdomain.

A further aspect of the invention comprises the minimization of a costfunctional using the approximation that the sum, weighted by windowfunctions, of the output of adjacent filters of the M filters is equalwhen applied to the same signal in time regions where said windowfunctions overlap. Preferably the method includes the step ofmultiplying M filtered estimates with temporal window functions. Theapplication of the temporal window functions, and hence the resultingtemporal windows, to the combined components ensures that the filteringprocess is local in time and allows the method adaptively to removenoise from the seismic data in accordance with a global optimizationcriterion. The data selection temporal window functions are preferablydetermined by two requirements, wherein the first requirement is thatthe sum over all windows at any given time equals unity and the secondrequirement is that only adjoining windows overlap. These requirementsensure that the global optimization of the filtered signal can be solvedby use of an approximation in which for the sum over all time and allfilters and all neighboring filters, the error function of a neighboringfilter is replaced with the error function associated with the filteritself.

The application of the data selection temporal windows decouples theequation required to solve the optimization of the filtered signal.

According to yet another aspect of the invention, the response of thefilter can be controlled by using a regularization parameter. Theparameter as applied herein determines the relative weight of twocomponents of the cost functional. One of the component of the costfunctional can be defined as output power, while the other can becharacterized as being essentially the white noise gain of the filterbank, i.e., the output of the filter in response to an inputuncorrelated in time and space.

The noisy signal may be pre-processed before being passed to theadaptive filtering means by dividing the signal into frequency bandsusing a reconstructing filter, for example a quadrature mirror filter.This allows a reduction in the number of data points to be processed andalso allows a reduction in the number of coefficients in the adaptivefiltering means as effectively reducing the bandwidth of the originalsignal.

The invention is applicable for two-dimensional (2D) andthree-dimensional (3D) seismic surveys, and can be used in land seismic,marine seismic including sea bottom seismic, and transitional zoneseismic.

The method can be performed on stored data or on raw seismic data as itis acquired. Thus raw seismic data may be filtered according to themethod at the data acquisition site. This ensures that a “cleaned”signal is available from the data acquisition site and may be downloadeddirectly from the site in this form. This reduces the amount of datasent for analysis off-site and reduces the costs and storage problemsassociated with accumulating sufficient quantities of noisy data foranalysis off-site. The method can be advantageously applied tosingle-sensor recordings, i.e. to recordings prior to any group formingwhich combines the signals of two or more seismic sensors.

Although the description of the present invention is based on seismicsignal processing, it can be applied to sonic signals as used forexample for well logging applications. Specific seismic applicationsinclude swell noise or streamer noise attenuation, including streamernoise attenuation in a cross-flow acquisition, attenuation of groundroll or mud roll or other coherent noise from marine, land, ortransition zone data, seismic interference canceling, i.e. filteringnoise using the full aperture of a multi-streamer array, which is eithertowed in the water or deployed at the sea-bottom , or removal ofsea-floor reflections from the notional source signature estimation, atechnique described for example in the European Patent ApplicationEP-A-066423. Other applications include noise suppression for variousborehole seismic exploration methods, where either noise or signal haspreferential directions. Known borehole seismic methods includeseismic-while-drilling (SWD), vertical seismic profiling (VSP),look-ahead and look-around while drilling.

These and other features of the invention, preferred embodiments andvariants thereof, possible applications and advantages will becomeappreciated and understood by those skilled in the art from thefollowing detailed description and drawings.

DRAWINGS

FIG. 1 illustrates general elements of a seismic land acquisition;

FIG. 2 shows a general block diagram of an adaptive beamformer inaccordance an example of the present invention;

FIG. 3 shows the example of a region preserved by an adaptive beamformerin accordance with an example of the present invention;

FIG. 4 shows steps of defining a constraint matrix in accordance with anexample of the present invention;

FIG. 5 illustrates steps of using the constaint matrix in a process offiltering seismic recordings in accordance with the present invention;

FIG. 6 shows the example of a preserved region including a noise regionsuppressed by an adaptive beamforming process in accordance with anexample of the present invention;

FIGS. 7A-C show examples of defining a preserved region or regions inaccordance with the present invention.

MODE(S) FOR CARRYING OUT THE INVENTION

In the following the basic concepts underlying the invention aredeveloped in detail.

A typical land seismic acquisition is illustrated in FIG. 1.

A source 10 is activated, thus generating seismic waves 11, i.e acousticwaves with frequencies of less than 500 Hz. The waves travel though theinterior of the Earth 12 and are reflected at various locations. Eventhough only one reflector 13 is shown, typically there are manyreflectors, each reflecting a fraction of the seismic wave back to thesurface. At the surface, seismic waves are recorded by seismic sensors14 (geophones). These sensors are spread along a line or in atwo-dimensional pattern.

As an example for a major source of noise, the travel path 15 of theso-called “ground roll” is shown. The ground roll is direct wave energywhich propagates within layers close to the surface. It is onedistinguishing feature of the ground roll to have different propagationcharacteristics than signals reflected from a deeper reflector layer:The ground roll reaches the sensors of the depicted sensor line oneafter the other. In contrast, the reflected seismic signal 16 from areflector at great depth reaches all the sensors 14 almostsimultaneously. When translated into the frequency-wavenumber domain,known as f−k domain, the desired seismic signal therefore usually lieswithin a narrow cone around the f-axis (equivalent to small values of k)whereas the ground roll tends to have larger values of k.

Referring now to FIG. 2, there is shown a general block diagram of anadaptive beamformer in accordance with the present invention. It isassumed the presence of K sensors located at r_(k) with k=1, . . . , K.Each sensor k records a signal g_(k)(n) with n=1, . . . , N. The letter‘n’ is used as an index on discrete time samples. The sampling intervalis Δt. The signals g_(k),(n) are beamsteered using delays τ_(k) towardsa general “signal direction”. This is the general direction from whichthe seismic signals are expected to arrive. The beamsteered channelsx_(k),(n) are processed by local multichannel adaptive filters toproduce the output signal: $\begin{matrix}{{y(n)} = {\sum\limits_{i = 1}^{M}\quad {\sum\limits_{k = 1}^{K}\quad {\sum\limits_{v = {- L_{1}}}^{L_{2}}\quad {{h_{i}(n)}w_{ikv}{x_{k}( {n - v} )}}}}}} & \lbrack 1\rbrack\end{matrix}$

where w_(ik)ν(t) are the adjustable coefficients of the adaptivefilters, h_(i)(n) are the windows applied at the output end, M is thenumber of local multichannel adaptive filters (or the number of outputwindows), and L=L₁+L₂+1 is the number of coefficients per channel. Hereand below, a bar under a letter denotes a vector (small letter) or amatrix (capital letter).

Equation [1] can be rewritten as a (windowed) sum over a scalar productusing a tap-input vector x(n) at time t defined as:

x(n)≡[x ₁(n+L ₁), . . . , x ₁(n−L ₂),

x ₂(n+L ₁), . . . , x ₂(n−L ₂),

x _(K) (n+L ₁), . . . , x _(K)(n−L ₂)]^(T)  [2]

and a tap-weight vector w_(i) defined as

 w _(i) ≡[w _(i1(−L) ₁ ₎ , . . . , w _(i1L) ₂ , w _(i2(−L) ₁ ₎ , . . . ,w _(i2L) ₂ , w _(iK(−L) ₁ ₎ , . . . , w _(iKL) ₂ ]^(T)  [3]

Using definitions [2] and [3], equation [1] becomes $\begin{matrix}{{y(n)} = {{\sum\limits_{i = 1}^{M}\quad {{h_{i}(n)}{\underset{\_}{w}}_{i}^{T}{{\underset{\_}{x}}_{k}(n)}}} = \quad {\sum\limits_{i = 1}^{M}\quad {{h_{i}(n)}{{\underset{\_}{x}}^{T}(n)}{{\underset{\_}{w}}_{i}.}}}}} & \lbrack 4\rbrack\end{matrix}$

Equations [1] and [4] describe how to find the beamformer or filter bankoutput once the M tap-weight vectors w_(i) have been specified. Thesevectors are computed as the solution of an optimization problem which isdescribed below.

The optimization problem is defined as $\begin{matrix}{{\min\limits_{{\underset{\_}{w}}_{1},\ldots \quad,{\underset{\_}{w}}_{M}}J} = {\min\limits_{{\underset{\_}{w}}_{1},\ldots \quad,{\underset{\_}{w}}_{M}}\{ {J_{1} + {\frac{\delta^{2}}{KL}J_{2}}} \}}} & \lbrack 5\rbrack\end{matrix}$

subject to constraints

C ^(T) w _(i) =f  [6]

where i=1,2, . . . , M and $\begin{matrix}{J_{1} = {\sum\limits_{n = 1}^{N}\quad {{y^{2}(n)}\quad {and}}}} & \lbrack 7\rbrack \\{{J_{2} = {\sum\limits_{i = 1}^{M}\quad {{{\underset{\_}{w}}_{i}}^{2}{\sum\limits_{n = 1}^{N}\quad {{h_{i}(n)}{{\underset{\_}{x}(n)}}^{2}}}}}},} & \lbrack 8\rbrack\end{matrix}$

KL is the total number of filter coefficients, and ||·|| denotes the L₂norm. This cost functional is a linear combination of the output powerof the beamformer (the first term in eq. [5]), and the so-called“white-noise gain” of the beamformer weighted by the input power (thesecond term in eq. [5]). The relative weight of the two terms isadjusted by the δ² term. Including the “white-noise gain” of thebeamformer in the cost functional is intended to increase the beamformerrobustness in the presence of signal modeling uncertainties (sometimescalled perturbations,) and numerical correlation between the signal andthe noise.

Equation [6] describes Q linear constraints on the admissible solutionsto the optimization problem. Here, the KLxQ matrix C is the constraintmatrix, and the Q-vector f is the response vector. The actual design ofthe linear constraints are discussed below.

A possible solution of the optimization depends on imposing thefollowing two constraints on the window functions h_(i)(n):$\begin{matrix}{{\sum\limits_{i = 1}^{M}\quad {h_{i}(n)}} = 1} & \lbrack 9\rbrack\end{matrix}$

for n =1, 2, . . . , N, and

h _(i)(n)h _(j)(n)=0  [10]

for j<>i−1,i,i+1. The first constraint ensures that the filter bank isequivalent to the single filter case if all the local filters (w_(i))areidentical. The second constraint ensures that the windows have compactsupport.

The optimization problem can be to a large extend decoupled using thesecond condition(eq. [10]), and the approximation $\begin{matrix}{{\sum\limits_{n}\quad {\sum\limits_{i}{\sum\limits_{{j = {i - 1}},{i + 1}}\quad {{h_{i}(n)}{h_{j}(n)}{\underset{\_}{w}}_{i}^{T}{\underset{\_}{x}(n)}{{\underset{\_}{x}}^{T}(n)}{\underset{\_}{w}}_{j}}}}} \approx {\sum\limits_{n}\quad {\sum\limits_{i}{\sum\limits_{{j = {i - 1}},{i + 1}}\quad {{h_{i}(n)}{h_{j}(n)}{\underset{\_}{w}}_{i}^{T}{\underset{\_}{x}(n)}{{\underset{\_}{x}}^{T}(n)}{{\underset{\_}{w}}_{i}.}}}}}} & \lbrack 11\rbrack\end{matrix}$

The approximation of equation [11] requires that neighboring filtersproduce similar results when applied to the same input data in timeregions where adjacent windows overlap, instead of requiring thatneighboring filters are similar on a point-by-point basis. Thus, theapproximation is similar to requiring that the integral of two functionsare close, rather than the functions themselves.

With this approximation, the first term of the cost functional, J_(i),becomes $\begin{matrix}{J_{1} = {\sum\limits_{i = 1}^{M}\quad {{\underset{\_}{w}}_{i}^{T}{\underset{\_}{\Phi}}_{i}{\underset{\_}{w}}_{i\quad}\quad {with}}}} & \lbrack 12\rbrack \\{{\underset{\_}{\Phi}}_{i} = {\sum\limits_{n}\quad {{h_{i}(n)}{\underset{\_}{x}(n)}{{{\underset{\_}{x}}^{T}(n)}.}}}} & \lbrack 13\rbrack\end{matrix}$

The second term in the cost functional can be rewritten as:$\begin{matrix}{{J_{2} = {\sum\limits_{i = 1}^{M}\quad {{w_{i}}^{2}{tr}\{ {\sum\limits_{n = 1}^{N}\quad {{h_{i}(n)}{\underset{\_}{x}(n)}{{\underset{\_}{x}}^{T}(n)}}} \}}}},} & \lbrack 14\rbrack\end{matrix}$

with “tr” denoting the trace of a matrix.

Combining Equations (5), (12), (14), and reorganizing the terms, thetotal cost functional can be written as: $\begin{matrix}{{J = {\sum\limits_{i = 1}^{M}\quad {{\underset{\_}{w}}_{i}^{T}\{ {{\underset{\_}{\Phi}}_{i} + {\frac{\delta^{2}}{KL}{{tr}( {\underset{\_}{\Phi}}_{i} )}\underset{\_}{I}}} \} {\underset{\_}{w}}_{i}}}},} & \lbrack 15\rbrack\end{matrix}$

where I denotes the KLxKL identity matrix. The decoupled optimizationproblem can be solved for each of the M time windows subject to theconstraints [6] Using the method of Lagrange multipliers, the optimaltap-weight in each window is given by

w _(i)*={tilde over (Φ)}_(i) ⁻¹ C(C ^(T){tilde over (Φ)}_(i) ⁻¹ C)⁻¹f,  [16]

with $\begin{matrix}{{{\underset{\_}{\overset{\sim}{\Phi}}}_{i} = {{\underset{\_}{\Phi}}_{i} + {\frac{\delta^{2}}{KL}\quad {{tr}( {\underset{\_}{\Phi}}_{i} )}\underset{\_}{I}}}},} & \lbrack 17\rbrack\end{matrix}$

The second term of the modified local correlation matrix Φ_(i) ⁻ can beregarded as a regularization term with δ² as the regularizationparameter. In array signal processing literature, regularization ofcorrelation matrices with the addition of a scaled identity matrix hasbeen suggested to increase robustness in the presence of perturbations,in the context of narrow-band beamforming. Here, the cost functional [5]includes the regularization term from the beginning leading to ageneralization for wide-band adaptive beamforming. Hence, the filterresponse changes as a function of the frequency of the signal.

When the input data to the beamformer is characterized by spatially andtemporally uncorrelated (or white) noise, both the correlation matrixΦ_(i) and the modified correlation matrix Φ_(i) ⁻ become proportional tothe identity matrix. In this case, the optimal weight vector becomes

w _(i) *≡w _(q) =C(C ^(T) C)⁻¹ f.  [18]

The weight vector w_(q) is called the quiescent solution to the optimalbeamformer problem, and the corresponding response is known as thequiescent response. Note that the quiescent solution depends entirely onthe constraint matrix C and the response vector f.

The optimal weight vector w_(i) approaches the quiescent weight vectorw_(q) even for general noise fields as the regularization parameter δ²is increased. In this case, the modified correlation matrix Φ_(i) ⁻approaches the identity matrix (cf. [17]). The regularization parameterδ² therefore weights the optimal solution between a solution that isentirely dependent on the received data, and a solution that isindependent of the data. For δ²=1, both solutions are equally weightedin the sense that their corresponding correlation matrices have equaltrace value. In situations where the perturbations are higher, i.e. theassumptions about the seismic acquisition geometry do not exactly hold,finding a beamformer response with a higher level of regularization cangive more robust results.

Another aspect of the invention relates to the design of linearconstraints (eq. [6]) to be imposed on the beamformer.

One type of linear constraints that can be imposed on the beamformer arethose designed to preserve seismic signals incident from a targetdirection, while suppressing interferences incident from otherdirections. Steering delays τ_(k) as those shown in FIG. 2 define asingle “look-direction”. Signals incident in this direction are inphase, and for these signals the system can be considered as a singleFIR (finite impulse response) filter. The values of the coefficients forthis equivalent processor are equal to the sums of the correspondingcoefficients in the adaptive processor. Each local beamformer w_(i)consists of the adaptive filters w_(i1), w_(i2), w_(iK) processing datafrom each channel, and a summing unit. The sum of the individual filtersw_(i1), w_(i2), . . . , w_(iK) is constrained to give w_(eq), which isthe desired response for signals incident in the look-direction, e.g., aunit pulse in the look direction. The quiescent response then becomesthat of a fixed-weight beamformer with single equal weights for allelements. In the frequency-wavenumber domain, this corresponds to a syncfunction that is constant in the f direction. Therefore, for increasingvalues of the regularization parameter δ², the beamformer preservessignals incident not only from the look direction, but also fromneighboring directions.

As discussed in the last section, using single look-directionconstraints and regularization, it is possible to preserve signalsincident from directions near the look direction. While this approach isuseful and sufficient for many applications, it is desirable to derivemore general linear constraints that will satisfy the requirements inany seismic data acquisition situation more directly.

In narrow-band beamforming, different generalized constraint designapproaches are known. Derivative constraints are used to influenceresponse over a region of the response space by forcing the derivativesof the beamformer response at certain points of the response space to bezero. Eigenvector constraints are based on a least squares approximationto the desired response, and are usually used to control the beamformerresponse over regions of the response space. Generalization of thesemethods to wide-band beamforming problems have shown that while theyprovide a good response in selected regions of the response space, theycan generate unacceptably high sidelobes in other regions.

For the present invention, the requirements of the generalizedconstraint design are to impose an arbitrary quiescent response on thebeamformer and to make sure that certain areas in thefrequency-wavenumber domain are entirely controlled by the quiescentresponse. These requirements have been established with the followingfunctional objectives in mind:

accommodation of an arbitrary range of apparent signal velocities;

increased robustness to perturbations;

capability to use larger arrays;

being able to run the adaptive beamformer with a lower regularizationlevel (δ²), hence achieving higher noise attenuation; and

achieving higher noise attenuation for a given level of regularizationby the appropriate design of the quiescent response.

To impose an arbitrary quiescent response on the beamformer, use can bemade of the fact that the linear constraints [6] define a Q-dimensionalhyperplane in a KL-dimensional space. Equation [18] shows that thequiescent weight vector w_(q) is the minimum norm solution to Equation[6], i.e., it is the shortest vector from the origin to the hyperplane.

Equation [18] also shows that w_(q) is a member of the subspace spannedby the columns of the constraint matrix C. The columns of C are ingeneral independent (otherwise some constraints would be redundant),thus they can be chosen to be orthogonal without loss of generality.After defining a desired quiescent weight vector w_(qd), this suggeststhe following forms for the constraint matrix C and the response vectorf:

C=[αw _(qd) , D]  [19]

and

f=[₀ ^(β)],  [20]

with the condition

β=α∥w _(qd)∥²  [21]

where D is a KL×(Q−1) matrix whose columns are orthogonal to w_(qd). Theexact form of the matrix D is described below. With C and f chosenaccording to [19] and [20], respectively, it can be shown that thedesired weight vector equals the quiescent response vector w_(q) (eq.[18]).

After defining the first column of the constraint matrix C and theresponse vector f to impose the quiescent weight vector, the definitionof the matrix D which is a part of C is derived in the followingsections.

In a seismic acquisition, reflection signals that should be preservedcan be considered as a linear combination of plane waves with associatedfrequency and wavenumber values from a known region of thefrequency-wavenumber space. This region, which is denoted A in FIG. 3,depends on the particular acquisition geometry, but is usually a conearound the frequency axis. One possible example of a preserved region inthe frequency-wavenumber domain is shown in FIG. 3, where A is chosen soas to include all signals of apparent velocity of +/−1500 m/s or more.In the present example, the beamformer response in region A should becontrolled entirely by a quiescent response which preserves the signal.

The set S_(A) of seismic signals to be preserved by the filteringprocess is defined by

S _(A) ={s(t, r): s(t, r)=∫∫_(A) dfdkS(f, k)e ^(j2π(ft−k·r))}  [22]

as composites of plane waves with associated frequencies and wavenumbervalues from a region A, where S(f,k) is the complex Fourier amplitudecorresponding to the plane wave component of a signal with frequency fand wavenumber k.

Using [22], the tap-input vector [2] can be rewritten as

s(n)=∫∫_((f,kεA)) dfdkS(f, k)e ^(j2πfnΔt) d(f, k),  [23]

with d(f,k) being defined as the array steering vector corresponding tothe plane wave component specified by particular frequency f andwavenumber k. It is noteworthy that in contrast to the example describedabove no time delays τ have been introduced into the signal path tosteer the filter response. The array steering vector can be written as aKronecker product: $\begin{matrix}{{\underset{\_}{d}( {f,\underset{\_}{k}} )} = {\begin{bmatrix}^{{- {j2\pi}}\quad {\underset{\_}{k} \cdot {\underset{\_}{r}}_{1}}} \\^{{- {j2\pi}}\quad {\underset{\_}{k} \cdot {\underset{\_}{r}}_{2}}} \\\vdots \\^{{- {j2\pi}}\quad {\underset{\_}{k} \cdot {\underset{\_}{r}}_{k}}}\end{bmatrix} \otimes {\begin{bmatrix}^{{j2\pi}\quad {fL}_{1}\Delta \quad t} \\\vdots \\1 \\\vdots \\^{{- {j2\pi}}\quad {fL}_{2}\Delta \quad t}\end{bmatrix}.}}} & \lbrack 24\rbrack\end{matrix}$

Using [4], the response of the beamformer to the signal tap-input vectors(n) is $\begin{matrix}{{y(n)} = {\sum\limits_{i = 1}^{M}\quad {{h_{i}(n)}{\int{\int_{A}{{f}\quad {\underset{\_}{k}}{S( {f,\underset{\_}{k}} )}^{T}^{{j2\pi}\quad {fn\Delta}\quad t}{{\underset{\_}{d}}^{T}( {f,\underset{\_}{k}} )}{{\underset{\_}{w}}_{i}.}}}}}}} & \lbrack 25\rbrack\end{matrix}$

For the beamformer response to be the same for both the optimal weightvector w_(i) and the quiescent weight vector w_(q), and furtherrequiring that this equality to hold for all signals s(t; r) of thepreserved region, i.e. signals with arbitrary associated Fouriercoefficients S(f, k) such that (f, k) is in A. This requires that

d ^(T)(f, k)w _(i) *=d ^(T)(f, k)w _(q), ∀(f, k)εA.  [26]

By decomposing the optimal weight vector into a fixed weight portionequal to the quiescent weight vector and an adaptive weight portionaccording to a solution known as “generalized sidelobe canceller” (GSC),it can be shown that the last equation is equivalent to requiring thatd(f,k) lies in the column space of the constraint matrix C.

It is therefore a further object of the present invention to find aefficient, i.e. preferably low rank, basis for the space of steeringvectors d(f,k). However, a scalar multiple of w_(qd) has already beeninstalled as the first column of C, we actually need to find a low rankbasis for the part of this space that lies in the orthogonal complementsubspace of w_(qd), The projection of d(f,k) onto the orthogonalcomplement of w_(qd) is the projected steering vector:

{tilde over (d)}(f, k)≡(I−P _(w) _(qd) )d(f, k),  [27]

where the expression in parentheses is the orthogonal complementprojector with respect to w_(qd) with

P _(w) _(qd) =w _(qd))w _(qd) ^(T) w _(qd))⁻¹ w _(qd) ^(T)  [28]

Using the fact that any KL-dimensional d⁻(f,k) can be written as alinear combination of orthonormal vectors {v_(1.), . . . , v_(XL)},$\begin{matrix}{{{\underset{\_}{\overset{\sim}{d}}( {f,\underset{\_}{k}} )} = {{\sum\limits_{p = 1}^{KL}\quad {{\sigma_{p}( {f,\underset{\_}{k}} )}{\underset{\_}{v}}_{p}}} \equiv {\underset{\_}{V\quad \sigma}( {f,k} )}}},} & \lbrack 29\rbrack\end{matrix}$

a rank P (P<KL ) approximation of the projected steering vectors isobtained by $\begin{matrix}{{{\overset{\sim}{{\underset{\_}{d}}_{p}}( {f,\underset{\_}{k}} )} = {{\sum\limits_{p = 1}^{P}\quad {{{\underset{\_}{\hat{\sigma}}}_{p}( {f,\underset{\_}{k}} )}{\underset{\_}{v}}_{p}}} \equiv {\underset{\_}{V}\quad {\hat{\underset{\_}{\sigma}}( {f,\underset{\_}{k}} )}}}},} & \lbrack 30\rbrack\end{matrix}$

where

 {circumflex over (σ)}(f, k)=[σ₁(f, k), . . . , σ_(p)(f, k), 0, . . . ,0]^(T).  [31]

To derive an efficient rank P representation of d⁻(f,k) for any (f,k) inregion A, an error functional with respect to the L₂ norm is defined as$\begin{matrix}{\mu_{P} \equiv {\int{\int_{A}\quad {{f}{\underset{\_}{k}}{{{{\underset{\_}{\overset{\sim}{d}}( {f,\underset{\_}{k}} )} - {{\underset{\_}{\overset{\sim}{d}}}_{P}( {f,\underset{\_}{k}} )}}}^{2}.}}}}} & \lbrack 32\rbrack\end{matrix}$

Using the correlation of all projected steering vectors in region A ofthe frequency-wavenumber space given by $\begin{matrix}{{\overset{\sim}{\underset{\_}{R}}}_{A} \equiv {\int{\int_{A}\quad {{f}{\underset{\_}{k}}{\underset{\_}{\overset{\sim}{d}}( {f,\underset{\_}{k}} )}{{\underset{\_}{\overset{\sim}{d}}}^{H}( {f,\underset{\_}{k}} )}}}}} & \lbrack 33\rbrack\end{matrix}$

the error functional can be expressed as $\begin{matrix}{\mu_{P} = {\sum\limits_{p = {P + 1}}^{KL}\quad {{\underset{\_}{v}}_{p}^{H}\overset{\sim}{\underset{\_}{R}}{{\underset{\_}{v}}_{p}.}}}} & \lbrack 34\rbrack\end{matrix}$

The superscript “H” denotes the conjugate transpose of a vector ormatrix.

The optimum set of ordered basis vectors v can be found by minimizingthe cost functional μ_(p) subject to the constraint that v_(p)^(H)v_(p)=1, with 1<=p<=KL. Using Lagrange multipliers, the task is tominimize $\begin{matrix}{\sum\limits_{p = {P + 1}}^{KL}{\lbrack \quad {{{\underset{\_}{v}}_{p}^{H}{\overset{\sim}{\underset{\_}{R}}}_{A}{\underset{\_}{v}}_{p}} - {\lambda_{p}( {{{\underset{\_}{v}}_{p}^{H}{\underset{\_}{v}}_{p}} - 1} )}} \rbrack.}} & \lbrack 35\rbrack\end{matrix}$

By taking the gradient with respect to v_(p) and setting it to zero, theoptimal basis vectors v₁, . . . , v_(KL) are found as the eigenvectorsof R_(A) ⁻(with respect to the eigenvalues λ_(p)). The missing part D ofthe constraints matrix C (cf. [19]) can now be defined as the principaleigenvectors of R_(A) ⁻:

D=[v ₁ , . . . , v _(p)].  [36]

Note that the steering vectors d(f; k) are in general complex valued.Therefore, their correlation matrix R_(A) ⁻ over a general region A inthe frequency-wavenumber space is complex valued, making theeigenvectors of R_(A) ⁻ hence the columns of C also complex valued.However, in seismics the signals are real valued signals which havecomplex conjugate Fourier coefficients. Therefore the types of A regionsthat are of interest are always symmetric in the frequency-wavenumberspace with respect to the origin. The resulting matrices (R_(A) ⁻, C)arethen all real valued.

The above described expansion of the projected steering vectors d⁻(f; k)is analogous to the Rarhunen-Lòeve expansion. While the originalKarhunen-Loeve expansion is for a random vector, the expansion presentedhere is for a deterministic set of vectors. This is reflected in the waythe approximation error functional μ_(p) is defined, cf. [32].

The covariance matrix of steering vectors, similar to the correlationmatrices defined in [33] was first introduced in by K. M. Buckley, IEEETrans. Acoust. Speech Signal Processing, Vol ASSP-35, 249-266, March1987, but was then heuristically defined within a stochastic framework,assuming zero mean signals and using a narrow-band representation ofwideband signals. In the description of the present invention, thecorrelation matrix has been derived from first principles within adeterministic framework.

The main steps of the generalized constraint design method are shown inthe flow diagram of FIG. 4. They include:

specification of the desired quiescent weight vector, w_(qd), whichdefines the first column of the constraint matrix;

specification of the signal protection region A in thefrequency-wavenumber space;

computation of R_(A) ⁻, the correlation matrix of all the projectedsteering vectors in region A; and

determination of the principal eigenvectors {v₁, . . . , v_(KL)} ofR_(A) ⁻ as the remaining columns of the constraint matrix.

Having computed these, the constraint matrix is specified as

C=└w _(qd) /∥w _(qd)∥² , v ₁ , . . . , v _(p)┘  [37]

and the response vector as $\begin{matrix}{{\underset{\_}{f} = \begin{bmatrix}1 \\\underset{\_}{0}\end{bmatrix}},} & \lbrack 38\rbrack\end{matrix}$

using a convenient choice for α and β in eq. [21]

The specification of the desired quiescent weight vector to form adesired quiescent response is essentially a non-adaptivemultidimensional filter design problem, for which many techniques exist.Reference can be made for example to handbooks such as W. Chen (ed.),“The Circuits and Filters Handbook”, IEEE and CRC Press, 2732-2761(1995), D. E. Dudgeon and R. M. Martinez, “Multidimensional DigitalSignal Processing”, Prentice Hall (1984), or J. S. Lim, Two-DimensionalSignal and Image Processing, Prentice Hall (1990).

Once the constraint matrix is defined, digitized recordings of singleseismic sensors can be filtered to generate a single sensor recordingwith reduced noise. The “cleaned” recording can then be used in furtherprocessing steps, such as group forming, stacking, velocity analysis,moveout correction etc., known in the art to ultimately generate arepresentation of subterranean formations. These steps are outlined inthe flow diagram of FIG. 5. As however the details of these steps (withthe exception of the filtering step) are of no particular concernregarding the present invention, are detailed description thereof isomitted herein.

The following section presents alternative methods of efficientlydefining the region A (see FIG. 3) protected by the quiescent response.

For seismic applications, the quiescent weight vector could be designedsuch that over the region A in the frequency-wavenumber space theresponse is near unity, thus preserving the seismic signals in thatregion. In regions of the frequency-wavenumber space where the noise isexpected to be present, the quiescent response should have low values,so that even when regularization is used, high performance can beachieved.

The constraint design process can be extended as described in the nextsection. The constraint design outlined above resulted in a low rankbasis of projected steering vectors in A. The objective was to preserveall signals in the preserved region A without any reference to theirrelative strength. This is reflected in the choice of the errorfunctional μ_(p), defined in [32]. In many applications this choicemakes sense, where it is desirable to protect signal components whichhave much lower amplitude then other signal components. On the otherhand, in some other applications it may be desirable to minimize thepower of the overall signal distortion.

This extension of the above described method is reflected by ageneralization of the definition of μ_(p) to

μ_(p)≡∫∫_(A) dfdkS(f, k)∥{tilde over (d)}(f, k)−{tilde over (d)} _(p)(f,k)∥².  [39]

The correlation matrix of the original steering vectors in A thenbecomes

R _(A)≡∫∫_(A) dfdk|S(f, k)|² d(f, k)d ^(H)(f, k).  [40]

A correlation matrix of the projected steering vectors can be derivedfrom [40] using the orthogonal complement projector (eqs. [27], [28])

{tilde over (R)} _(A)=(I−P _(w) _(qd) )R _(A)(I−P _(w) _(qd) ).  [41]

In some applications, it may be desirable to reserve additional regionsof the frequency-wavenumber space over which the beamformer responsewould be controlled entirely by the quiescent response. In FIG. 6, thereis shown an example of region A in the frequency-wavenumber space whichincludes signal protection region A1 and noise region A2. If it is knownthat in a certain environment there is almost always a coherent noisecomponent occupying the region A2, then it may be beneficial to designthe quiescent weight vector w_(qd) which would put deep nulls in thatregion, i.e, which suppresses any signal from a direction (f,k) where(f,k) is element of A2. Then, the adaptive weights would concentrate onattenuating the noise in the remaining regions of thefrequency-wavenumber space.

Employing the generalized constraint design methodology, the quiescentresponse of the beamformer can be specified exactly, and the beamformerresponse over a pre-specified region A in the f−k space can beconstrained to approximate the quiescent response to an arbitraryextent. The accuracy of this approximation is controlled by the userparameter P, which is the number of principal components of R_(A) ⁻. AsP is increased, more and more degrees of freedom of the beamformer arefixed, and the adaptive degrees of freedom are reduced. As P approachesKL−1 the beamformer response approaches the quiescent response not onlyin A, but over the entire f−k space regardless of the noise fieldcharacteristics or the regularization parameter used.

When the region A is relatively large, the distribution of theeigenvalues of R_(A) ⁻ may be only slowly decreasing, so that the numberof principal components required to adequately represent the space ofsteering vectors in A may be high. If only a small number of principalcomponents are used in order to keep more degrees of freedom adaptive,the response of the beamformer may deviate from the quiescent responsein A significantly. This may happen even over the original lookdirection indicated by k=0. This is in contrast with the original single“look-direction” constraints described before, which guarantee that theresponse over the k=0 line is identically unity for all frequencies.

In order to guarantee that the response of the beamformer is exactlythat of the quiescent response over k=0 as in the case of singlelook-direction constraints while constraining the response arbitrarilytightly in the rest of region A, the region A is preferably partitionedand each section is treated separately.

Methods of pursuing this approach are shown in FIGS. 7A-C, where A isthe region as shown in FIG. 3. In the following examples, this region Ais partioned into sections A1 and A2 as illustrated in FIGS. 7A and 7B.

It can be shown that if A1 is the k=0 line in the f−k space, then thespace of projected steering vectors over A1 can be spanned by L−1eigenvectors of R_(A1) ⁻(see eq. [42] below), as long as the quiescentresponse w_(g) has unity response over A1, as usually is the case. Thatis, the matrix $\begin{matrix}{{\overset{\sim}{\underset{\_}{R}}}_{A1} \equiv {\int_{A1}\quad {{f}{\underset{\_}{\overset{\sim}{d}}( {f,\underset{\_}{0}} )}{{{\underset{\_}{\overset{\sim}{d}}}^{H}( {f,\underset{\_}{0}} )}.}}}} & \lbrack 42\rbrack\end{matrix}$

has rank L−1. This is an interesting observation, as it neatly ties thesingle look-direction constraints and the generalized constraints. Ifthe signal protection region A is chosen as A1, the quiescent responsew_(q) is chosen as $\begin{matrix}{{{\underset{\_}{w}}_{q} = {\begin{bmatrix}1 \\1 \\\vdots \\1\end{bmatrix} \otimes \begin{bmatrix}0 \\\vdots \\\frac{1}{K} \\\vdots \\0\end{bmatrix}}},} & \lbrack 43\rbrack\end{matrix}$

and the columns of matrix D in the constraint matrix C are chosen as theL−1 eigenvectors of R_(A1) ⁻ corresponding to the non-zero eigenvalues,then the linear constraints set by the two methods become identical.Having above defined the region A1 as the k axis, there are variouspossibilities for defining A2. One would be to define A2 as the rest ofthe region A as in FIG. 7A. Another way to define A2 is shown in FIG.7B. In this case A2 is basically the boundary of the original region A.

In general, if the number of array elements K is small and the originalregion A has a relatively short wavenumber extent, setting theadditional constraints as in FIG. 7B would be preferred. When the numberof sensor elements is small, the variability of the beamformer responseas a function of wavenumber is limited, and constraining the responsemore strictly at a few points for each frequency may be more effectivethan constraining it over all the points of a wide region. If the numberof array elements is large or the region A is wide in wavenumbers,defining A2 as in FIG. 7A would be preferred. It is of course possibleto generalize even further to combine the two types of additionalconstraints as shown in FIG. 7C giving rise to three sections denotedA1, A2 and A3 of region A.

When the protection region A is partitioned using any of the approachesabove, the constraints can be computed defining the response vector asin eq. [38], where 0 has appropriate length, and the constraint matrix Cis given as

C=└w _(qd) /∥w _(qd)∥² , D _(A1) , D _(A2)┘  [44]

Here, D_(A1) is the matrix whose columns are the L−1 principaleigenvectors of R_(A1) ⁻ (cf. eq.[42])

{tilde over (R)} _(A1)=(I−P _(w) _(qd) )R _(A1)(I−P _(w) _(qd) )  [45]

using the projection matrix with respect to w_(qd) defined in eqs. [27],[28]. D_(A2) is the matrix whose columns are the principal eigenvectorsof R_(A2) ⁻:

{tilde over (R)} _(A2)=(I−P _(C) _(A1) )R _(A2)(I−P _(C) _(A1) )=(I−P_(D) _(A1) )(I−P _(w) _(qd) )R _(A2)(I−P _(w) _(qd) )(I−P _(D) _(A1))′  [46]

where

C _(A1) =└w _(qd) /∥w _(qd)∥² ,D _(A1)┘,  [47]

P _(C) _(A1) =C _(A1)(C _(A1) ^(T) C _(A1))⁻¹ C _(A1) ^(T)  [48]

P _(D) _(A1) =D _(A1)(D _(A1) ^(T) D _(A1))⁻¹ D _(A1) ^(T), and  [49]

R _(A2)=∫∫_(A2) dfdkd(f, k)d ^(H)(f, k).  [50]

Equation [146] shows two alternative ways of computing R_(A1) ⁻. Theseequations have been written for partitioning of A into two subregions,but can be generalized to more sub-regions.

By partitioning the original signal protection region A as described,choosing A1 as the k=0 line, setting the quiescent response w_(q) as in[43], and constructing the columns of matrix D in the constraint matrixC are chosen as the L−1 eigenvectors of R_(A1) ⁻ corresponding to thenonzero eigenvalues, any additional sub-regions A2, A3, etc. would giveadditional control over beamformer response with respect to thelook-direction constraints.

For some applications, it may be useful to reduce the degrees of freedomused by the adaptive beamformer. In the so-called partially adaptivebeamformer, only a portion of the available degrees of freedom are usedadaptively. The main advantages of reducing the adaptive degrees offreedom are reduced computational cost and improved adaptive convergencerate. The primary disadvantage of partially adaptive beamforming is adegradation in the steady state interference cancellation capability ofthe beamformer. Therefore, the objective of partially adaptivebeamformer design is to reduce the number of adaptive weights withoutsignificantly degrading the performance of the adaptive array.

Previous partially adaptive methods includes numerical techniques forapproximately minimizing the average generalized sidelobe canceller(GSC) output power for a desired number of adaptive weights, where theaverage is over a range of jammer parameters. The present invention usesa method which is based on a design method described by H. Yang and M.A. Ingram, IEEE Trans. On Antennas and Propagation, Vol. 45, 843-850,May 1997. It also attempts to minimize the average GSC output power, butunder a constraint that the reduced-dimensional solutions for all of thescenario trials lie in the same subspace. This constraint makes itpossible to use a singular value decomposition to get the rank-reducingtransformation, thereby simplifying the optimization problem.

The generalized sidelobe canceller solution can be written as (cf.[18]):

w _(i) *=w _(q) −Bw _(ai),  [51]

where B is a KL×(KL−Q) full rank matrix whose columns span theorthogonal subspace of the constraints matrix C and is known as theblocking matrix. The vector w_(ai) is the KL×Q dimensional adaptive partof the optimal weight vector and is given by

w _(ai)=(B ^(T){tilde over (Φ)}_(i) B)⁻¹ B ^(T){tilde over (Φ)}_(i) w_(q).  [52]

The partially adaptive GSC achieves a smaller number W of adaptiveweights, through the use of a (KL−Q)×W linear transformation T followingB prior to adaptive weighting. The partially adaptive optimal weightvector can be expressed as

w _(i) *=w _(q) −BTw _(pi),  [53]

where the W-dimensional adaptive part of the optimal weight is

w _(pi)=(T ^(T) B ^(T){tilde over (Φ)}_(i) BT)⁻¹ T ^(T) B ^(T){tildeover (Φ)}_(i) w _(q).  [54]

It is now the aim to choose T which minimizes the interference and noiseoutput power over a set of likely interference scenarios. Thesescenarios can be characterized by different parameters like the numberof interferers, interferer directions, interferer spectral densities,white noise levels, etc. The applied method can be summarized asfollows:

for each random outcome θ_(j) from a distribution of scenarioparameters, compute the full-rank optimal adaptive weight vector w_(ai)from [52] and the transformed weight vector α given by

α(θ_(j))=UΣU ^(T) w _(ai)(θ_(j)),  [55]

where

B ^(T){tilde over (Φ)}_(i)(θ_(j))B=UΣ ² U ^(T)  [56]

is the eigendecomposition of Φ_(i) ⁻(θj);

store vectors w_(ai)(θ_(j))and α(θ_(j)) into the matrices W and A,respectively;

compute the singular value decomposition of A to get U_(A) from

A=U _(A)Σ_(A) V _(A) ^(H);  [57]

and

derive T as the first W columns of WA^(#)U_(A), where the superscript“#” indicates the pseudoinverse.

In most seismic surveys, the noise such as ground roll or swell noiseoccupies only a fraction of the temporal bandwidth available. Forexample in a land seismic survey, the Nyquist frequency is 250 Hz, whilemost of the ground roll energy is under 30 Hz. Concentrating filteringefforts to the frequency band where the noise resides is desirable toreduce computational cost.

One means of achieving this aim involves adding QMF (quadrature mirrorfilter) perfect reconstruction filter banks, described for example by P.P. Vaidyanathan, in “Multirate Systems and Filter Banks, Prentice Hall,1993 or by M. J. T. Smith and T. P. Barnwell III, in: IEEE Trans.Acoust. Speech Signal Processing, Vol. ASSP-34, 434-441 (1986) to theseismic noise and interference suppression system using adaptivemultichannel filter banks, as described above. Two filter banks are usedin this system. The QMF filter bank is used to decompose the traces intofrequency bands and decimate before adaptive filtering is applied, andis subsequently used for resynthesizing the original signal. Themultichannel adaptive filter bank is the heart of the system performingthe actual filtering for noise suppression. Using the perfectreconstruction filter banks to decimate reduces the number of points tobe processed and also allows reduction in the number of coefficients inthe adaptive filters, bringing in significant savings in CPU time andmemory requirements.

What is claimed is:
 1. Method for filtering noise from discrete noisyseismic signals, said method comprising the steps of: receiving signalsrepresenting seismic energy reflected and/or refracted from an earthstructure using a plurality of seismic receivers; and filtering receivedsignals using an at least partially adaptive filter that attenuatessignals that lack predetermined propagation characteristics; saidfiltering step comprising the steps of: defining at least twoindependent sets of conditions with a first set defining a desiredresponse, and a second set, defining the propagation characteristics ofsignals to be preserved wherein the first set is not affected by thesecond set; and adapting filter coefficients of said filter subject tosaid at least two independent sets of conditions so as to optimize thefilter output for signals that lack the propagation characteristics ofsignals to be preserved.
 2. The method of claim 1 wherein the step ofadapting comprises optimizing a cost function representing the outputpower of the filter and the white noise gain of the filter.
 3. Themethod of claim 2 wherein the cost function represents a weighted sum ofthe output power of the filter and the white noise gain of the filter.4. The method of claim 3 wherein the relative weight of the output powerand the “white noise gain” in the cost function is adjustable.
 5. Themethod of claim 1 wherein the propagation characteristics of a signal isdetermined by its temporal and spatial spectral content.
 6. The methodof claim 1 wherein the propagation characteristics of a signal isdetermined by its wavenumber and frequency.
 7. The method of claim 1wherein the propagation characteristics of a signal is defined as aregion or as regions in the frequency-wavenumber domain.
 8. The methodof claim 1 wherein the propagation characteristics of a signal includesa plurality directions in the frequency-wavenumber domain.
 9. The methodof claim 1 wherein the conditions define a constraint matrix comprisingat least two mutually orthogonal subspaces with a first subspacedefining the desired response of the filter and a second subspacedefining the propagation characteristics of the signals to be preservedwith said matrix being applied during the adaptation process of filtercoefficients.
 10. The method of claim 1 wherein the conditions define aconstraint matrix using a quiescent response vector and principleeigenvectors of a correlation matrix made of steering vectors definingat least one region of the frequency-wavenumber domain wherein thesignals are preserved.
 11. The method of claim 1 wherein the conditionsare defined such that signals with pre-determined propagationcharacteristics are suppressed by the filter.
 12. The method of claim 1,wherein said filtering step further comprises defining further sets ofconditions with one set of conditions forcing the desired response ofthe filter to zero and one set of conditions defining the propagationcharacteristics of signals to be suppressed by the filter.
 13. Themethod of claim 1 wherein the filter comprises M temporally localfilters, said filters forming a filter bank, and M being a number equalto or larger than two.
 14. The method of claim 13 including the step ofmultiplying M filtered estimates with temporal window functions.
 15. Themethod of claim 14 wherein the temporal window functions arecharacterized by the requirement that only adjoining windows overlap.16. The method of claim 14 wherein the cost function is minimized usingthe approximation that the sum, weighted by window functions, of theoutput of adjacent filters of the M filters is equal when applied to thesame signal in time regions where said window functions overlap.
 17. Themethod of claim 1 comprising the further step of dividing filtercoefficients into a fixed part and a adaptive part.
 18. The method ofclaim 1 comprising the further step of splitting the signals into atleast two frequency bands, filtering at least one of said at least twofrequency band, and, before further processing, recombining saidfrequency bands to regain noise attenuated seismic signals.
 19. Themethod of claim 1 wherein the discrete noisy signals used as input arerecordings from individual seismic sensors prior to any group formingtechniques.
 20. The method of claim 1 comprising the further steps ofpositioning at least one seismic source and a plurality of seismicreceivers at appropriate positions in a land, marine or transitionalzone environment; activating said at least one source to propagateenergy through subterranean formations, and using said receivers tomeasure energy as single sensor recordings, converting said singlesensor recordings into discrete seismic signals; and transmitting saidsignals as input to the filter.
 21. A method for filtering noise fromdiscrete noisy seismic signals, said method comprising the steps of:receiving signals representing seismic energy reflected and/or refractedfrom an earth structure using a plurality of seismic receivers; andfiltering received signals using an at least partially adaptive filterthat attenuates signals that lack certain propagation characteristics;said filtering step comprising the steps of: defining conditions with afirst set defining a desired response, and a second set definingpropagation characteristics of signals to be preserved; and adaptingfilter coefficients of said filter subject to said sets of conditions soas to optimize the filter output for signals that lack the propagationcharacteristics of signals to be preserved at least in party byoptimizing a cost function representing the output power of the filterand the white noise gain of the filter.
 22. The method of claim 21wherein the cost function represents a weighted sum of the output powerof the filter and the white noise gain of the filter.
 23. The method ofclaim 22 wherein the relative weight of the output power and the “whitenoise gain” in the cost function is adjustable.